import xarray as xr
import numpy as np
import netCDF4 as nc
from scipy import signal
import numpy.ma as ma
import regionmask
from scipy import stats
import scipy
import regionmask
from scipy.interpolate import interp2d
from scipy import interpolate
from cartopy.util import add_cyclic_point
import matplotlib.pyplot as plt
import matplotlib as mpl
from matplotlib.colors import ListedColormap, LinearSegmentedColormap
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import cartopy.mpl.ticker as cticker
import cartopy.io.shapereader as shpreader

treefrac2 = xr.open_dataset('D:/Python/XXHH/gpp-ssp126/gpp_Lmon_CanESM5_ssp126_r1i1p2f1_gn_201501-210012.nc')
treefrac22 = xr.open_dataset('D:/Python/XXHH/gpp-ssp126-ssp370/gpp_Lmon_CanESM5_ssp126-ssp370Lu_r1i1p2f1_gn_201501-210012.nc')
treefrac1 = xr.open_dataset('D:/Python/XXHH/gpp-ssp126/gpp_Lmon_ACCESS-ESM1-5_ssp126_r1i1p1f1_gn_201501-210012.nc')
treefrac11 = xr.open_dataset('D:/Python/XXHH/gpp-ssp126-ssp370/gpp_Lmon_ACCESS-ESM1-5_ssp126-ssp370Lu_r1i1p1f1_gn_201501-210012.nc')
treefrac3 = xr.open_dataset('D:/Python/XXHH/gpp-ssp126/gpp_Lmon_CMCC-ESM2_ssp126_r1i1p1f1_gn_201501-210012.nc')
treefrac33 = xr.open_dataset('D:/Python/XXHH/gpp-ssp126-ssp370/gpp_Lmon_CMCC-ESM2_ssp126-ssp370Lu_r1i1p1f1_gn_201501-210012.nc')
treefrac4 = xr.open_dataset('D:/Python/XXHH/gpp-ssp126/gpp_Lmon_IPSL-CM6A-LR_ssp126_r1i1p1f1_gr_201501-210012.nc')
treefrac44 = xr.open_dataset('D:/Python/XXHH/gpp-ssp126-ssp370/gpp_Lmon_IPSL-CM6A-LR_ssp126-ssp370Lu_r1i1p1f1_gr_201501-210012.nc')
treefrac5 = xr.open_dataset('D:/Python/XXHH/gpp-ssp126/gpp_Lmon_MPI-ESM1-2-LR_ssp126_r1i1p1f1_gn_201501-210012.nc')
treefrac55 = xr.open_dataset('D:/Python/XXHH/gpp-ssp126-ssp370/gpp_Lmon_MPI-ESM1-2-LR_ssp126-ssp370Lu_r1i1p1f1_gn_201501-210012.nc')
treefrac6 = xr.open_dataset('D:/Python/XXHH/gpp-ssp126/gpp_Lmon_GFDL-ESM4_ssp126_r1i1p1f1_gr1_201501-210012.nc')
treefrac66 = xr.open_dataset('D:/Python/XXHH/gpp-ssp126-ssp370/gpp_Lmon_GFDL-ESM4_ssp126-ssp370Lu_r1i1p1f1_gr1_201501-210012.nc')
print(treefrac6)
treeFrac1 = treefrac1.gpp.loc['2015-01-01':'2100-12-16', :, :]
treeFrac2 = treefrac2.gpp.loc['2015-01-01':'2100-12-16', :, :]
treeFrac3 = treefrac3.gpp.loc['2015-01-01':'2100-12-16', :, :]
treeFrac4 = treefrac4.gpp.loc['2015-01-01':'2100-12-16', :, :]
treeFrac5 = treefrac5['gpp']
treeFrac6 = treefrac6['gpp']
treeFrac11 = treefrac11.gpp.loc['2015-01-01':'2100-12-16', :, :]
treeFrac22 = treefrac22.gpp.loc['2015-01-01':'2100-12-16', :, :]
treeFrac33 = treefrac33.gpp.loc['2015-01-01':'2100-12-16', :, :]
treeFrac44 = treefrac44.gpp.loc['2015-01-01':'2100-12-16', :, :]
treeFrac55 = treefrac55['gpp']
treeFrac66 = treefrac66['gpp']
print(treeFrac55)

lon1 = treeFrac1.lon

lat1 = treeFrac1.lat
print(lat1)
print(lon1)
lon2 = treeFrac2.lon
lat2 = treeFrac2.lat
lon3 = treeFrac3.lon
lat3 = treeFrac3.lat
lon4 = treeFrac4.lon
lat4 = treeFrac4.lat
lon5 = treeFrac5.lon
lat5 = treeFrac5.lat
lon6 = treeFrac6.lon
lat6 = treeFrac6.lat

treeFrac1 = np.array(treeFrac1).reshape((86,12,145,192))
treeFrac2 = np.array(treeFrac2).reshape((86,12,64,128))
treeFrac3 = np.array(treeFrac3).reshape((86,12,192,288))
treeFrac4 = np.array(treeFrac4).reshape((86,12,143,144))
treeFrac5 = np.array(treeFrac5).reshape((86,12,96,192))
treeFrac6 = np.array(treeFrac6).reshape((86,12,180, 288))

treeFrac11 = np.array(treeFrac11).reshape((86,12,145,192))
treeFrac22 = np.array(treeFrac22).reshape((86,12,64,128))
treeFrac33 = np.array(treeFrac33).reshape((86,12,192,288))
treeFrac44 = np.array(treeFrac44).reshape((86,12,143,144))
treeFrac55 = np.array(treeFrac55).reshape((86,12,96,192))
treeFrac66 = np.array(treeFrac66).reshape((86,12,180, 288))

treeFrac1 = np.nanmean(treeFrac1, axis=1)
print(treeFrac1.shape)
treeFrac2 = np.nanmean(treeFrac2, axis=1)
treeFrac3 = np.nanmean(treeFrac3, axis=1)
treeFrac4 = np.nanmean(treeFrac4, axis=1)
treeFrac5 = np.nanmean(treeFrac5, axis=1)
treeFrac6 = np.nanmean(treeFrac6, axis=1)
treeFrac11 = np.nanmean(treeFrac11, axis=1)
print(treeFrac1.shape)
treeFrac22 = np.nanmean(treeFrac22, axis=1)
treeFrac33 = np.nanmean(treeFrac33, axis=1)
treeFrac44 = np.nanmean(treeFrac44, axis=1)
treeFrac55 = np.nanmean(treeFrac55, axis=1)
treeFrac66 = np.nanmean(treeFrac66, axis=1)

tree11 = treeFrac11 - treeFrac1
tree22 = treeFrac22 - treeFrac2
tree33 = treeFrac33 - treeFrac3
tree44 = treeFrac44 - treeFrac4
tree55 = treeFrac55 - treeFrac5
tree66 = treeFrac66 - treeFrac6

tree11 = tree11*1000*60*60*24*365
tree22 = tree22*1000*60*60*24*365
tree33 = tree33*1000*60*60*24*365
tree44 = tree44*1000*60*60*24*365
tree55 = tree55*1000*60*60*24*365
tree66 = tree66*1000*60*60*24*365


Q1 = np.load('D:/Python/XXHH/area/ACCESS-ESM1-5-area.npz')
Q2 = np.load('D:/Python/XXHH/area/CanESM5-area.npz')
Q3 = np.load('D:/Python/XXHH/area/CMCC-ESM2-area.npz')
Q4 = np.load('D:/Python/XXHH/area/IPSL-CM6A-LR-area.npz')
Q5 = np.load('D:/Python/XXHH/area/MPI-ESM1-2-LR-area.npz')
Q6 = np.load('D:/Python/XXHH/area/GFDL-ESM4-area.npz')

area1 = Q1['area']
area2 = Q2['area']
area3 = Q3['area']
area4 = Q4['area']
area5 = Q5['area']
area6 = Q6['area']

tree11 = tree11 * area1 / 10**9
tree22 = tree22 * area2 / 10**9
tree33 = tree33 * area3 / 10**9
tree44 = tree44 * area4 / 10**9
tree55 = tree55 * area5 / 10**9
tree66 = tree66 * area6 / 10**9


Q1 = np.load('D:/Python/XXHH/area/ACCESS-ESM1-5-forest.npz')
Q2 = np.load('D:/Python/XXHH/area/CanESM5-forest.npz')
Q3 = np.load('D:/Python/XXHH/area/CMCC-ESM2-forest.npz')
Q4 = np.load('D:/Python/XXHH/area/IPSL-CM6A-LR-forest.npz')
Q5 = np.load('D:/Python/XXHH/area/MPI-ESM1-2-LR-forest.npz')
Q6 = np.load('D:/Python/XXHH/area/GFDL-ESM4-forest.npz')

forest1 = Q1['quyu']
forest2 = Q2['quyu']
forest3 = Q3['quyu']
forest4 = Q4['quyu']
forest5 = Q5['quyu']
forest6 = Q6['quyu']

sorted_list1 = np.zeros((86,145,192))
#min_values1 = np.zeros((9,145,192))

sorted_list2 = np.zeros((86,64,128))

sorted_list3 = np.zeros((86,192,288))

sorted_list4 = np.zeros((86,143,144))

sorted_list5 = np.zeros((86,96,192))

sorted_list6 = np.zeros((86,180,288))

for j in range(0, 145):
    for k in range(0, 192):
        sorted_list1[:,j,k]= sorted(forest1[:,j,k])
        s1 = max(1, int(len(sorted_list1[:,j,k]) * 0.1))


for j in range(0, 64):
    for k in range(0, 128):
        sorted_list2[:,j,k]= sorted(forest2[:,j,k])
        s2 = max(1, int(len(sorted_list2[:,j,k]) * 0.1))

for j in range(0, 192):
    for k in range(0, 288):
        sorted_list3[:,j,k]= sorted(forest3[:,j,k])
        s3 = max(1, int(len(sorted_list3[:, j, k]) * 0.1))

for j in range(0, 143):
    for k in range(0, 144):
        sorted_list4[:,j,k]= sorted(forest4[:,j,k])
        s4 = max(1, int(len(sorted_list4[:, j, k]) * 0.1))

for j in range(0, 96):
    for k in range(0, 192):
        sorted_list5[:,j,k]= sorted(forest5[:,j,k])
        s5 = max(1, int(len(sorted_list5[:, j, k]) * 0.1))


for j in range(0, 180):
    for k in range(0, 288):
        sorted_list6[:,j,k]= sorted(forest6[:,j,k])
        s6 = max(1, int(len(sorted_list6[:, j, k]) * 0.1))


index1 = np.zeros((s1,145,192))
GPP1 = np.zeros((s1,145,192))
forest_index1 = np.zeros((s1,145,192))
gpp_forest1 = np.zeros((145,192))

index2 = np.zeros((s2,64,128))
GPP2 = np.zeros((s2,64,128))
forest_index2 = np.zeros((s2,64,128))
gpp_forest2 = np.zeros((64,128))

index3 = np.zeros((s3,192,288))
GPP3 = np.zeros((s3,192,288))
forest_index3 = np.zeros((s3,192,288))
gpp_forest3 = np.zeros((192,288))

index4 = np.zeros((s4,143,144))
GPP4 = np.zeros((s4,143,144))
forest_index4 = np.zeros((s4,143,144))
gpp_forest4 = np.zeros((143,144))

index5 = np.zeros((s5,96,192))
GPP5 = np.zeros((s5,96,192))
forest_index5 = np.zeros((s5,96,192))
gpp_forest5 = np.zeros((96,192))

index6 = np.zeros((s6,180,288))
GPP6 = np.zeros((s6,180,288))
forest_index6 = np.zeros((s6,180,288))
gpp_forest6 = np.zeros((180,288))


for j in range(0, 145):
    for k in range(0, 192):
        #index1[:,j,k]= np.argsort(tree11[:,j,k])
        #GPP[:,j,k] = index1[:9,j,k]
        index1[:,j,k] = np.argpartition(tree11[:,j,k] , s1)[:s1]
        GPP1[:,j,k] = tree11[[np.argpartition(tree11[:,j,k] , s1)[:s1], j, k]]
        forest_index1[:,j,k] = np.argpartition(forest1[:,j,k] , s1)[:s1]
print(index1.shape)
print(GPP1.shape)
print(forest_index1.shape)

for j in range(0, 64):
    for k in range(0, 128):
        index2[:,j,k] = np.argpartition(tree22[:,j,k] , s2)[:s2]
        GPP2[:,j,k] = tree22[[np.argpartition(tree22[:,j,k] , s2)[:s2], j, k]]
        forest_index2[:,j,k] = np.argpartition(forest2[:,j,k] , s2)[:s2]

for j in range(0, 192):
    for k in range(0, 288):
        index3[:,j,k] = np.argpartition(tree33[:,j,k] , s3)[:s3]
        GPP3[:,j,k] = tree33[[np.argpartition(tree33[:,j,k] , s3)[:s3], j, k]]
        forest_index3[:,j,k] = np.argpartition(forest3[:,j,k] , s3)[:s3]

for j in range(0, 143):
    for k in range(0, 144):
        index4[:,j,k] = np.argpartition(tree44[:,j,k] , s4)[:s4]
        GPP4[:,j,k] = tree44[[np.argpartition(tree44[:,j,k] , s4)[:s4], j, k]]
        forest_index4[:,j,k] = np.argpartition(forest4[:,j,k] , s4)[:s4]

for j in range(0, 96):
    for k in range(0, 192):
        index5[:,j,k] = np.argpartition(tree55[:,j,k] , s5)[:s5]
        GPP5[:,j,k] = tree55[[np.argpartition(tree55[:,j,k] , s5)[:s5], j, k]]
        forest_index5[:,j,k] = np.argpartition(forest5[:,j,k] , s5)[:s5]

for j in range(0, 180):
    for k in range(0, 288):
        index6[:,j,k] = np.argpartition(tree66[:,j,k] , s6)[:s6]
        GPP6[:,j,k] = tree66[[np.argpartition(tree66[:,j,k] , s6)[:s6], j, k]]
        forest_index6[:,j,k] = np.argpartition(forest6[:,j,k] , s6)[:s6]

for j in range(0, 145):
    for k in range(0, 192):
        for i in range(0, 8):
            if index1[i, j, k] in forest_index1[:, j, k]:
                gpp_forest1[j, k] = gpp_forest1[j, k] + GPP1[i, j, k]
            else:
                continue
print(gpp_forest1.shape)

for j in range(0, 64):
    for k in range(0, 128):
        for i in range(0, 8):
            if index2[i, j, k] in forest_index2[:, j, k]:
                gpp_forest2[j, k] = gpp_forest2[j, k] + GPP2[i, j, k]
            else:
                continue

for j in range(0, 192):
    for k in range(0, 288):
        for i in range(0, 8):
            if index3[i, j, k] in forest_index3[:, j, k]:
                gpp_forest3[j, k] = gpp_forest3[j, k] + GPP3[i, j, k]
            else:
                continue

for j in range(0, 143):
    for k in range(0, 144):
        for i in range(0, 8):
            if index4[i, j, k] in forest_index4[:, j, k]:
                gpp_forest4[j, k] = gpp_forest4[j, k] + GPP4[i, j, k]
            else:
                continue

for j in range(0, 96):
    for k in range(0, 192):
        for i in range(0, 8):
            if index5[i, j, k] in forest_index5[:, j, k]:
                gpp_forest5[j, k] = gpp_forest5[j, k] + GPP5[i, j, k]
            else:
                continue

for j in range(0, 180):
    for k in range(0, 288):
        for i in range(0, 8):
            if index6[i, j, k] in forest_index6[:, j, k]:
                gpp_forest6[j, k] = gpp_forest6[j, k] + GPP6[i, j, k]
            else:
                continue



Q1 = np.load('D:/Python/XXHH/area/ACCESS-ESM1-5-tree-50.npz')
Q2 = np.load('D:/Python/XXHH/area/CanESM5-tree-50.npz')
Q3 = np.load('D:/Python/XXHH/area/CMCC-ESM2-tree-50.npz')
Q4 = np.load('D:/Python/XXHH/area/IPSL-CM6A-LR-tree-50.npz')
Q5 = np.load('D:/Python/XXHH/area/MPI-ESM1-2-LR-tree-50.npz')
Q6 = np.load('D:/Python/XXHH/area/GFDL-ESM4-tree-50.npz')

quyu1 = Q1['quyu']
quyu2 = Q2['quyu']
quyu3 = Q3['quyu']
quyu4 = Q4['quyu']
quyu5 = Q5['quyu']
quyu6 = Q6['quyu']


gpp_forest1 = np.where(quyu1<0, gpp_forest1, np.nan)
gpp_forest1 = np.where(np.isnan(gpp_forest1), 0, gpp_forest1)
gpp_forest2 = np.where(quyu2<0, gpp_forest2, np.nan)
gpp_forest2 = np.where(np.isnan(gpp_forest2), 0, gpp_forest2)
gpp_forest3 = np.where(quyu3<0, gpp_forest3, np.nan)
gpp_forest3 = np.where(np.isnan(gpp_forest3), 0, gpp_forest3)
gpp_forest4 = np.where(quyu4<0, gpp_forest4, np.nan)
gpp_forest4 = np.where(np.isnan(gpp_forest4), 0, gpp_forest4)
gpp_forest5 = np.where(quyu5<0, gpp_forest5, np.nan)
gpp_forest5 = np.where(np.isnan(gpp_forest5), 0, gpp_forest5)
gpp_forest6 = np.where(quyu6<0, gpp_forest6, np.nan)
gpp_forest6 = np.where(np.isnan(gpp_forest6), 0, gpp_forest6)

lon1 = np.array(lon1)
lat1 = np.array(lat1)
lon2 = np.array(lon2)
lat2 = np.array(lat2)
lon3 = np.array(lon3)
lat3 = np.array(lat3)
lon4 = np.array(lon4)
lat4 = np.array(lat4)
lon5 = np.array(lon5)
lat5 = np.array(lat5)
lon6 = np.array(lon6)
lat6 = np.array(lat6)
lon_new = np.arange(0,360,1)
lat_new = np.arange(-90,90,1)

f1 = interpolate.interp2d(lon1, lat1, gpp_forest1, kind='linear')
min_values1_new = f1(lon_new, lat_new)
print(min_values1_new.shape)

f2 = interp2d(lon2, lat2, gpp_forest2, kind='linear')
min_values2_new = f2(lon_new, lat_new)
print(min_values2_new)

f3 = interp2d(lon3, lat3, gpp_forest3, kind='linear')
min_values3_new = f3(lon_new, lat_new)
print(min_values3_new)

f4 = interp2d(lon4, lat4, gpp_forest4, kind='linear')
min_values4_new = f4(lon_new, lat_new)
print(min_values4_new)

f5 = interp2d(lon5, lat5, gpp_forest5, kind='linear')
min_values5_new = f5(lon_new, lat_new)
print(min_values5_new)

f6 = interp2d(lon6, lat6, gpp_forest6, kind='linear')
min_values6_new = f6(lon_new, lat_new)
print(min_values6_new)


min_values1_new = np.where(min_values1_new==0, np.nan, min_values1_new)
min_values2_new = np.where(min_values2_new==0, np.nan, min_values2_new)
min_values3_new = np.where(min_values3_new==0, np.nan, min_values3_new)
min_values4_new = np.where(min_values4_new==0, np.nan, min_values4_new)
min_values5_new = np.where(min_values5_new==0, np.nan, min_values5_new)
min_values6_new = np.where(min_values6_new==0, np.nan, min_values6_new)
#min_values1 = np.where(quyu1<0, min_values1, np.nan)

min_values1_new = min_values1_new*1000
min_values2_new = min_values2_new*1000
min_values3_new = min_values3_new*1000
min_values4_new = min_values4_new*1000
min_values5_new = min_values5_new*1000
min_values6_new = min_values6_new*1000

cycle_sst1, cycle_lon1 = add_cyclic_point(min_values1_new, coord=lon_new)
cycle_LON1, cycle_LAT1 = np.meshgrid(cycle_lon1, lat_new)
cycle_sst2, cycle_lon2 = add_cyclic_point(min_values2_new, coord=lon_new)
cycle_LON2, cycle_LAT2 = np.meshgrid(cycle_lon2, lat_new)
cycle_sst3, cycle_lon3 = add_cyclic_point(min_values3_new, coord=lon_new)
cycle_LON3, cycle_LAT3 = np.meshgrid(cycle_lon3, lat_new)
cycle_sst4, cycle_lon4 = add_cyclic_point(min_values4_new, coord=lon_new)
cycle_LON4, cycle_LAT4 = np.meshgrid(cycle_lon4, lat_new)
cycle_sst5, cycle_lon5 = add_cyclic_point(min_values5_new, coord=lon_new)
cycle_LON5, cycle_LAT5 = np.meshgrid(cycle_lon5, lat_new)
cycle_sst6, cycle_lon6 = add_cyclic_point(min_values6_new, coord=lon_new)
cycle_LON6, cycle_LAT6 = np.meshgrid(cycle_lon6, lat_new)

mean_modle = np.zeros((6,180,361))
mean_modle[0, :, :] = cycle_sst1
mean_modle[1, :, :] = cycle_sst2
mean_modle[2, :, :] = cycle_sst3
mean_modle[3, :, :] = cycle_sst4
mean_modle[4, :, :] = cycle_sst5
mean_modle[5, :, :] = cycle_sst6

mean_modle = np.nanmean(mean_modle, axis=0)

f_g = nc.Dataset('D:/Python/XXHH/deforestion_extreme_GPP_mean_modle.nc','w',format = 'NETCDF4')


f_g.createDimension('lat',len(lat_new))
f_g.createDimension('lon',len(cycle_lon1))


##创建变量。参数依次为：‘变量名称’，‘数据类型’，‘基础维度信息’

f_g.createVariable('lat',np.float32,('lat'))
f_g.createVariable('lon',np.float32,('lon'))
f_g.createVariable('gpp', np.float32, ('lat','lon'))
#f_g.createVariable('GPP', np.float32, ('lat','lon'))


lon = np.array(cycle_lon1)
print(lon.shape)
f_g.variables['lon'][:] = lon
lat = np.array(lat_new)
f_g.variables['lat'][:] = lat
var_data = mean_modle
f_g.variables['gpp'][:, :] = var_data

f_LR1 = xr.open_dataset('D:/Python/XXHH/deforestion_extreme_GPP_mean_modle.nc')
print(f_LR1)

Q1 = np.load('D:/Python/XXHH/area/GPP-mean_modle.npz')
GPP_mean_modle = Q1['mean_modle']

#percent = (GPP_mean_modle - mean_modle) / GPP_mean_modle
#percent = (GPP_mean_modle - mean_modle)
percent = mean_modle / GPP_mean_modle
percent = percent*100
percent1 = percent
percent1 = np.nanmean(percent1,axis=(0,1))
print(percent1)

np.savez('D:/Python/XXHH/area/Fig_2b.npz', Region=mean_modle)
np.savez('D:/Python/XXHH/area/Fig_2d.npz', Region=percent)


proj = ccrs.PlateCarree(central_longitude=0)
#leftlon, rightlon, lowerlat, upperlat = (-180,180,-90,90)
leftlon, rightlon, lowerlat, upperlat = (0,361,-90,90)
img_extent = [leftlon, rightlon, lowerlat, upperlat]
fig = plt.figure(figsize=(12,8))

#ax1 = fig.add_axes([0.05, 0.74, 0.22, 0.32],projection = proj) 一列四行
ax1 = fig.add_axes([0.05, 0.7, 0.22, 0.32],projection = proj)
#ax1.set_title('(a)',loc='left',fontsize=22,family='Times New Roman')
ax1.set_title('(a) ACCESS-ESM1-5',loc='left',fontsize=11,family='Times New Roman')
ax1.set_extent(img_extent, crs=ccrs.PlateCarree())
ax1.add_feature(cfeature.COASTLINE.with_scale('50m'))
ax1.add_feature(cfeature.LAKES, alpha=0.5)
#ax1.set_xticks(np.array([-180,-135,-90,-45,0,45,90,135,180]), crs=ccrs.PlateCarree())
ax1.set_xticks(np.array([-180,-120,-60,0,60,120,180]), crs=ccrs.PlateCarree())
ax1.set_yticks(np.array([-90,-60,-30,0,30,60,90]), crs=ccrs.PlateCarree())
ax1.xaxis.set_major_formatter(cticker.LongitudeFormatter())
ax1.yaxis.set_major_formatter(cticker.LatitudeFormatter())
plt.xticks(fontsize=9,family='Times New Roman')
plt.yticks(fontsize=9,family='Times New Roman')
cmap = mpl.cm.YlOrRd_r
#colors = ['#9B4739','#CE5D45','#E89767','#FBDA99','#E0EFAC','#A0D279','#76A459']
#cmap1 = LinearSegmentedColormap.from_list("mycmap", colors)
newcolors=cmap(np.linspace(0, 1, 256))
newcmp = ListedColormap(newcolors[0:256])
ax1.add_feature(cfeature.NaturalEarthFeature('physical', 'ocean', '110m', edgecolor='black', facecolor='white'))
#levels = (0,50,100,125,150,175, 200,225,250,275,300, 600)
#c1 = ax.contourf(X, Y, BPRE, zorder=0, extend = 'both',levels=levels, transform=ccrs.PlateCarree(), cmap=newcmp)
c1 = ax1.contourf(cycle_LON1, cycle_LAT1, cycle_sst1 , zorder=0, extend = 'both',transform=ccrs.PlateCarree(),levels=np.arange(-38,0,2),cmap=newcmp)
#ax = srex.plot(projection=ccrs.PlateCarree(central_longitude=0), add_label=False)
#ax1.add_geometries(china, ccrs.PlateCarree(), facecolor='none', edgecolor='k', linewidth=1, zorder=1)
#cb = plt.colorbar(c1, orientation='horizontal',fraction=0.05)

#ax2 = fig.add_axes([0.05, 0.47, 0.22, 0.32],projection = proj)
ax2 = fig.add_axes([0.05, 0.40, 0.22, 0.32],projection = proj)
#ax1.set_title('(a)',loc='left',fontsize=22,family='Times New Roman')
ax2.set_title('(c) CanESM5',loc='left',fontsize=11,family='Times New Roman')
ax2.set_extent(img_extent, crs=ccrs.PlateCarree())
ax2.add_feature(cfeature.COASTLINE.with_scale('50m'))
ax2.add_feature(cfeature.LAKES, alpha=0.5)
#ax1.set_xticks(np.array([-180,-135,-90,-45,0,45,90,135,180]), crs=ccrs.PlateCarree())
ax2.set_xticks(np.array([-180,-120,-60,0,60,120,180]), crs=ccrs.PlateCarree())
ax2.set_yticks(np.array([-90,-60,-30,0,30,60,90]), crs=ccrs.PlateCarree())
ax2.xaxis.set_major_formatter(cticker.LongitudeFormatter())
ax2.yaxis.set_major_formatter(cticker.LatitudeFormatter())
plt.xticks(fontsize=9,family='Times New Roman')
plt.yticks(fontsize=9,family='Times New Roman')
cmap = mpl.cm.YlOrRd_r
#colors = ['#9B4739','#CE5D45','#E89767','#FBDA99','#E0EFAC','#A0D279','#76A459']
#cmap1 = LinearSegmentedColormap.from_list("mycmap", colors)
newcolors=cmap(np.linspace(0, 1, 256))
newcmp = ListedColormap(newcolors[0:256])
ax2.add_feature(cfeature.NaturalEarthFeature('physical', 'ocean', '110m', edgecolor='black', facecolor='white'))
#levels = (0,50,100,125,150,175, 200,225,250,275,300, 600)
#c1 = ax.contourf(X, Y, BPRE, zorder=0, extend = 'both',levels=levels, transform=ccrs.PlateCarree(), cmap=newcmp)
c2 = ax2.contourf(cycle_LON2, cycle_LAT2, cycle_sst2 , zorder=0, extend = 'both',transform=ccrs.PlateCarree(),levels=np.arange(-38,0,2),cmap=newcmp)
#ax1.add_geometries(china, ccrs.PlateCarree(), facecolor='none', edgecolor='k', linewidth=1, zorder=1)
#cb = plt.colorbar(c2, orientation='horizontal',fraction=0.05)

#ax3 = fig.add_axes([0.05, 0.20, 0.22, 0.32],projection = proj)
ax3 = fig.add_axes([0.05, 0.10, 0.22, 0.32],projection = proj)
#ax1.set_title('(a)',loc='left',fontsize=22,family='Times New Roman')
ax3.set_title('(e) CMCC-ESM2',loc='left',fontsize=11,family='Times New Roman')
ax3.set_extent(img_extent, crs=ccrs.PlateCarree())
ax3.add_feature(cfeature.COASTLINE.with_scale('50m'))
ax3.add_feature(cfeature.LAKES, alpha=0.5)
#ax1.set_xticks(np.array([-180,-135,-90,-45,0,45,90,135,180]), crs=ccrs.PlateCarree())
ax3.set_xticks(np.array([-180,-120,-60,0,60,120,180]), crs=ccrs.PlateCarree())
ax3.set_yticks(np.array([-90,-60,-30,0,30,60,90]), crs=ccrs.PlateCarree())
ax3.xaxis.set_major_formatter(cticker.LongitudeFormatter())
ax3.yaxis.set_major_formatter(cticker.LatitudeFormatter())
plt.xticks(fontsize=9,family='Times New Roman')
plt.yticks(fontsize=9,family='Times New Roman')
cmap = mpl.cm.YlOrRd_r
#colors = ['#9B4739','#CE5D45','#E89767','#FBDA99','#E0EFAC','#A0D279','#76A459']
#cmap1 = LinearSegmentedColormap.from_list("mycmap", colors)
newcolors=cmap(np.linspace(0, 1, 256))
newcmp = ListedColormap(newcolors[0:256])
ax3.add_feature(cfeature.NaturalEarthFeature('physical', 'ocean', '110m', edgecolor='black', facecolor='white'))
#levels = (0,50,100,125,150,175, 200,225,250,275,300, 600)
#c1 = ax.contourf(X, Y, BPRE, zorder=0, extend = 'both',levels=levels, transform=ccrs.PlateCarree(), cmap=newcmp)
c3 = ax3.contourf(cycle_LON3, cycle_LAT3, cycle_sst3 , zorder=0, extend = 'both',transform=ccrs.PlateCarree(),levels=np.arange(-38,0,2),cmap=newcmp)
#ax1.add_geometries(china, ccrs.PlateCarree(), facecolor='none', edgecolor='k', linewidth=1, zorder=1)
#cb = plt.colorbar(c3, orientation='horizontal',fraction=0.05)


#ax4 = fig.add_axes([0.05, -0.07, 0.22, 0.32],projection = proj)
ax4 = fig.add_axes([0.32, 0.70, 0.22, 0.32],projection = proj)
#ax1.set_title('(a)',loc='left',fontsize=22,family='Times New Roman')
ax4.set_title('(b) IPSL-CM6A-LR',loc='left',fontsize=11,family='Times New Roman')
ax4.set_extent(img_extent, crs=ccrs.PlateCarree())
ax4.add_feature(cfeature.COASTLINE.with_scale('50m'))
ax4.add_feature(cfeature.LAKES, alpha=0.5)
#ax1.set_xticks(np.array([-180,-135,-90,-45,0,45,90,135,180]), crs=ccrs.PlateCarree())
ax4.set_xticks(np.array([-180,-120,-60,0,60,120,180]), crs=ccrs.PlateCarree())
ax4.set_yticks(np.array([-90,-60,-30,0,30,60,90]), crs=ccrs.PlateCarree())
ax4.xaxis.set_major_formatter(cticker.LongitudeFormatter())
ax4.yaxis.set_major_formatter(cticker.LatitudeFormatter())
plt.xticks(fontsize=9,family='Times New Roman')
plt.yticks(fontsize=9,family='Times New Roman')
cmap = mpl.cm.YlOrRd_r
#colors = ['#9B4739','#CE5D45','#E89767','#FBDA99','#E0EFAC','#A0D279','#76A459']
#cmap1 = LinearSegmentedColormap.from_list("mycmap", colors)
newcolors=cmap(np.linspace(0, 1, 256))
newcmp = ListedColormap(newcolors[0:256])
ax4.add_feature(cfeature.NaturalEarthFeature('physical', 'ocean', '110m', edgecolor='black', facecolor='white'))
#levels = (0,50,100,125,150,175, 200,225,250,275,300, 600)
#c1 = ax.contourf(X, Y, BPRE, zorder=0, extend = 'both',levels=levels, transform=ccrs.PlateCarree(), cmap=newcmp)
c4 = ax4.contourf(cycle_LON4, cycle_LAT4, cycle_sst4 , zorder=0, extend = 'both',transform=ccrs.PlateCarree(),levels=np.arange(-38,0,2),cmap=newcmp)
#ax1.add_geometries(china, ccrs.PlateCarree(), facecolor='none', edgecolor='k', linewidth=1, zorder=1)
#cb = plt.colorbar(c4, orientation='horizontal',fraction=0.05)


ax5 = fig.add_axes([0.32, 0.40, 0.22, 0.32],projection = proj)
#ax1.set_title('(a)',loc='left',fontsize=22,family='Times New Roman')
ax5.set_title('(d) MPI-ESM1-2-LR',loc='left',fontsize=11,family='Times New Roman')
ax5.set_extent(img_extent, crs=ccrs.PlateCarree())
ax5.add_feature(cfeature.COASTLINE.with_scale('50m'))
ax5.add_feature(cfeature.LAKES, alpha=0.5)
#ax1.set_xticks(np.array([-180,-135,-90,-45,0,45,90,135,180]), crs=ccrs.PlateCarree())
ax5.set_xticks(np.array([-180,-120,-60,0,60,120,180]), crs=ccrs.PlateCarree())
ax5.set_yticks(np.array([-90,-60,-30,0,30,60,90]), crs=ccrs.PlateCarree())
ax5.xaxis.set_major_formatter(cticker.LongitudeFormatter())
ax5.yaxis.set_major_formatter(cticker.LatitudeFormatter())
plt.xticks(fontsize=9,family='Times New Roman')
plt.yticks(fontsize=9,family='Times New Roman')
cmap = mpl.cm.YlOrRd_r
#colors = ['#9B4739','#CE5D45','#E89767','#FBDA99','#E0EFAC','#A0D279','#76A459']
#cmap1 = LinearSegmentedColormap.from_list("mycmap", colors)
newcolors=cmap(np.linspace(0, 1, 256))
newcmp = ListedColormap(newcolors[0:256])
ax5.add_feature(cfeature.NaturalEarthFeature('physical', 'ocean', '110m', edgecolor='black', facecolor='white'))
#levels = (0,50,100,125,150,175, 200,225,250,275,300, 600)
#c1 = ax.contourf(X, Y, BPRE, zorder=0, extend = 'both',levels=levels, transform=ccrs.PlateCarree(), cmap=newcmp)
c5 = ax5.contourf(cycle_LON5, cycle_LAT5, cycle_sst5 , zorder=0, extend = 'both',transform=ccrs.PlateCarree(),levels=np.arange(-38,0,2),cmap=newcmp)
#ax1.add_geometries(china, ccrs.PlateCarree(), facecolor='none', edgecolor='k', linewidth=1, zorder=1)
#cb = plt.colorbar(c5, orientation='horizontal',fraction=0.05)


ax6 = fig.add_axes([0.32, 0.10, 0.22, 0.32],projection = proj)
#ax1.set_title('(a)',loc='left',fontsize=22,family='Times New Roman')
ax6.set_title('(f) GFDL-ESM4',loc='left',fontsize=11,family='Times New Roman')
ax6.set_extent(img_extent, crs=ccrs.PlateCarree())
ax6.add_feature(cfeature.COASTLINE.with_scale('50m'))
ax6.add_feature(cfeature.LAKES, alpha=0.5)
#ax1.set_xticks(np.array([-180,-135,-90,-45,0,45,90,135,180]), crs=ccrs.PlateCarree())
ax6.set_xticks(np.array([-180,-120,-60,0,60,120,180]), crs=ccrs.PlateCarree())
ax6.set_yticks(np.array([-90,-60,-30,0,30,60,90]), crs=ccrs.PlateCarree())
ax6.xaxis.set_major_formatter(cticker.LongitudeFormatter())
ax6.yaxis.set_major_formatter(cticker.LatitudeFormatter())
plt.xticks(fontsize=9,family='Times New Roman')
plt.yticks(fontsize=9,family='Times New Roman')
cmap = mpl.cm.YlOrRd_r
#colors = ['#9B4739','#CE5D45','#E89767','#FBDA99','#E0EFAC','#A0D279','#76A459']
#cmap1 = LinearSegmentedColormap.from_list("mycmap", colors)
newcolors=cmap(np.linspace(0, 1, 256))
newcmp = ListedColormap(newcolors[0:256])
ax6.add_feature(cfeature.NaturalEarthFeature('physical', 'ocean', '110m', edgecolor='black', facecolor='white'))
#levels = (0,50,100,125,150,175, 200,225,250,275,300, 600)
#c1 = ax.contourf(X, Y, BPRE, zorder=0, extend = 'both',levels=levels, transform=ccrs.PlateCarree(), cmap=newcmp)
c6 = ax6.contourf(cycle_LON6, cycle_LAT6, cycle_sst6 , zorder=0, extend = 'both',transform=ccrs.PlateCarree(),levels=np.arange(-38,0,2),cmap=newcmp)
#ax1.add_geometries(china, ccrs.PlateCarree(), facecolor='none', edgecolor='k', linewidth=1, zorder=1)
#cb = plt.colorbar(c6, orientation='horizontal',fraction=0.05)
#cb.set_label('Tree Fraction(%)', fontsize=12,family='Times New Roman')
cb = fig.colorbar(c6, ax=[ax1, ax2, ax3, ax4, ax5, ax6],  shrink=0.6, orientation='horizontal',pad=0.07)
cb.set_label('$\Delta$GPP(Tg C)', fontsize=10,family='Times New Roman')
plt.savefig('D:/picture/second/figs8.pdf',dpi=400,bbox_inches='tight')
plt.show()

proj = ccrs.PlateCarree(central_longitude=0)
#leftlon, rightlon, lowerlat, upperlat = (-180,180,-90,90)
leftlon, rightlon, lowerlat, upperlat = (0,361,-90,90)
img_extent = [leftlon, rightlon, lowerlat, upperlat]
fig1 = plt.figure(figsize=(12,8))

#ax1 = fig.add_axes([0.05, 0.74, 0.22, 0.32],projection = proj) 一列四行
ax1 = fig1.add_axes([0.1, 0.5, 0.32, 0.42],projection = proj)
#ax1.set_title('(a)',loc='left',fontsize=22,family='Times New Roman')
ax1.set_title('(b) Negative extreme GPP under deforestation',loc='left',fontsize=15,family='Times New Roman')
ax1.set_extent(img_extent, crs=ccrs.PlateCarree())
ax1.add_feature(cfeature.COASTLINE.with_scale('50m'))
ax1.add_feature(cfeature.LAKES, alpha=0.5)
#ax1.set_xticks(np.array([-180,-135,-90,-45,0,45,90,135,180]), crs=ccrs.PlateCarree())
ax1.set_xticks(np.array([-180,-120,-60,0,60,120,180]), crs=ccrs.PlateCarree())
ax1.set_yticks(np.array([-90,-60,-30,0,30,60,90]), crs=ccrs.PlateCarree())
ax1.xaxis.set_major_formatter(cticker.LongitudeFormatter())
ax1.yaxis.set_major_formatter(cticker.LatitudeFormatter())
plt.xticks(fontsize=12,family='Times New Roman')
plt.yticks(fontsize=12,family='Times New Roman')
cmap = mpl.cm.YlOrRd_r
#colors = ['#9B4739','#CE5D45','#E89767','#FBDA99','#E0EFAC','#A0D279','#76A459']
#cmap1 = LinearSegmentedColormap.from_list("mycmap", colors)
newcolors=cmap(np.linspace(0, 1, 256))
newcmp = ListedColormap(newcolors[0:256])
ax1.add_feature(cfeature.NaturalEarthFeature('physical', 'ocean', '110m', edgecolor='black', facecolor='white'))
#levels = (0,50,100,125,150,175, 200,225,250,275,300, 600)
#c1 = ax.contourf(X, Y, BPRE, zorder=0, extend = 'both',levels=levels, transform=ccrs.PlateCarree(), cmap=newcmp)
c1 = ax1.contourf(cycle_LON1, cycle_LAT1, mean_modle, zorder=0, extend = 'both',transform=ccrs.PlateCarree(),levels=np.arange(-38,0,2),cmap=newcmp)
#ax = srex.plot(projection=ccrs.PlateCarree(central_longitude=0), add_label=False)
#ax1.add_geometries(china, ccrs.PlateCarree(), facecolor='none', edgecolor='k', linewidth=1, zorder=1)
cb = plt.colorbar(c1, orientation='horizontal',fraction=0.1)
cb.set_label('$\Delta$GPP(Tg C)', fontsize=12,family='Times New Roman')
#plt.savefig('D:/picture/fig 10b.tif',dpi=400,bbox_inches='tight')
plt.show()

proj = ccrs.PlateCarree(central_longitude=0)
#leftlon, rightlon, lowerlat, upperlat = (-180,180,-90,90)
leftlon, rightlon, lowerlat, upperlat = (0,361,-90,90)
img_extent = [leftlon, rightlon, lowerlat, upperlat]
fig1 = plt.figure(figsize=(12,8))

#ax1 = fig.add_axes([0.05, 0.74, 0.22, 0.32],projection = proj) 一列四行
ax1 = fig1.add_axes([0.1, 0.5, 0.32, 0.42],projection = proj)
#ax1.set_title('(a)',loc='left',fontsize=22,family='Times New Roman')
ax1.set_title('(d) Changes in GPP under deforestation',loc='left',fontsize=15,family='Times New Roman')
ax1.set_extent(img_extent, crs=ccrs.PlateCarree())
ax1.add_feature(cfeature.COASTLINE.with_scale('50m'))
ax1.add_feature(cfeature.LAKES, alpha=0.5)
#ax1.set_xticks(np.array([-180,-135,-90,-45,0,45,90,135,180]), crs=ccrs.PlateCarree())
ax1.set_xticks(np.array([-180,-120,-60,0,60,120,180]), crs=ccrs.PlateCarree())
ax1.set_yticks(np.array([-90,-60,-30,0,30,60,90]), crs=ccrs.PlateCarree())
ax1.xaxis.set_major_formatter(cticker.LongitudeFormatter())
ax1.yaxis.set_major_formatter(cticker.LatitudeFormatter())
plt.xticks(fontsize=12,family='Times New Roman')
plt.yticks(fontsize=12,family='Times New Roman')
cmap = mpl.cm.YlGn
#clist=['orange','honeydew','lightgreen','limegreen','mediumseagreen','green','forestgreen','seagreen', 'darkgreen']
# N 表示插值后的颜色数目
#newcmp = LinearSegmentedColormap.from_list('chaos',clist,N=256)
newcolors=cmap(np.linspace(0, 1, 256))
newcmp = ListedColormap(newcolors[0:200])
ax1.add_feature(cfeature.NaturalEarthFeature('physical', 'ocean', '110m', edgecolor='black', facecolor='white'))
#levels = (0,50,100,125,150,175, 200,225,250,275,300, 600)
#c1 = ax.contourf(X, Y, BPRE, zorder=0, extend = 'both',levels=levels, transform=ccrs.PlateCarree(), cmap=newcmp)
c1 = ax1.contourf(cycle_LON1, cycle_LAT1, percent, zorder=0, extend = 'both',transform=ccrs.PlateCarree(),levels=np.arange(-10,50,1),cmap=newcmp)
#ax = srex.plot(projection=ccrs.PlateCarree(central_longitude=0), add_label=False)
#ax1.add_geometries(china, ccrs.PlateCarree(), facecolor='none', edgecolor='k', linewidth=1, zorder=1)
cb = plt.colorbar(c1, orientation='horizontal',fraction=0.1)
cb.set_label('%', fontsize=12,family='Times New Roman')
#plt.savefig('D:/picture/fig 10d-1.tif',dpi=400,bbox_inches='tight')
plt.show()

